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THE FINITE ELEMENT METHOD IN 
LOW SPEED AERODYNAMICS 

by 

A. J. Baker 
Visiting Professor* 

Old Dominion University 

and 

Paul D. Manhardt 
Textron Bell Aerospace 


SUMMARY 

Within recent years, the finite element method has emerged as an 
alternative theoretical foundation for establishing numerical solution 
algorithms for field problems in continuum mechanics, and in particular 
fluid dynamics. The fundamental concept of the finite element pro- 
cedure is identification of a computational control -volume within which 
the conservation laws of mechanics are numerically approximated using 
a well-defined and uniformly consistent procedure. In contrast to the 
more familiar finite difference technique, the finite element algorithm 
explicitly retains use of the calculus and vector field theory in forma- 
tion and evaluation of the discretized-equi valent matrix statement of 
control -volume conservation properties. For specifically non-linear 
differential equations, like the Navier-Stokes system in fluid dynamics, 


*0n leave from Textron Bell Aerospace, 1974-75 academic year. 


the finite element method always yields a unique discretized-equivalent 
matrix expression. Interestingly, the resultant expressions appear uni- 
formly comparable to the computationally preferred forms established via 
trial and evaluation techniques using difference algebra. These theoret- 
ical features, coupled with its intrinsic capability to employ speci- 
fiably non-regular computational meshes, and to enforce non-trivial grad- 
ient-type boundary conditions on non-coordinate solution domain boundaries, 
appears to render the finite element procedure of particular potential use 
in computational aerodynamics. 

The results of a study, conducted for verification of this hypothesis 
and reported herein, indeed show that the finite element procedure can 
be of significant impact in design of the "computational wind tunnel" for 
low speed aerodynamics. The uniformity of the mathematical differential 
equation description, for viscous and/or inviscid, multi-dimensional sub- 
sonic flows about practical aerodynamic system configurations, is utilized 
to establish the general form of the finite element algorithm, universally 
applicable to all problem classes. The COMOC computer program, under 
development for several years as a finite element test bed, is similarly 
applicable to each of these diverse classes. Following completion of the 
theoretical developments, example numerical results for inviscid flow 
analyses, as well as viscous boundary-layer, parabolic- and full- Navier- 
Stokes flow descriptions, verify the capabilities and overall versatility 
of the fundamental algorithm for aerodynamics. The proven mathematical 
basis, coupled with the distinct user-orientation features of the computer 
program embodiment, portend near-term evolution of a highly useful analyt- 
ical design tool to support computational configuration studies in low 
speed aerodynamics. 


INTRODUCTION 


Development of the techniques that has emerged as 11 the finite element 
method, 11 began in the civil and structural engineering community in the early 
1950’s. Conventional engineering structures can be visualized as comprised of 
discrete elements interconnected at a finite number of points. For example, 
the analysis of pin-connected structures via imposition of a global force bal- 
ance., can determine tension and compression members and loads; it is a well 
established study at the undergraduate level. In an elastic (or any otherl) 
continuum, however, the number of such connections becomes infinite, and 
therein lies the difficulty in establishing a tractible engineering analysis. 
Turner et al, (ref, 1), in their original concept of a "finite element," 
attempted to bridge this gap by introducing a method to transform a continuum 
into an equivalent discretized finite assemblage of nodal behavior. The 
method developed rapidly thereafter, in good part due to the parallel develop- 
ment of large capacity digital computers. Existence of this computer hardware 
made possible the embodiment of the analytical developments into a viable tool 
(computer program) for engineering analysis of complex non-conventional 
structures. 

While a detailed static force balance formed the theoretical foundation 
for the earlier finite element developments, this was rapidly replaced by the 
energy concepts that also form a significant branch of conventional structural 
analysis. From this viewpoint, the "finite element method" emerged as a tech- 
nique for modeling the strain energy of a continuum structural system in terms 
of the behavior of local discrete subsystems. Total system energy is approxi- 
mately determined by the assembly of the incremental work done by the surface 
tractions and imbedded stress distribution as the equivalent discretized struc 
ture is loaded. The actual node point displacement distribution is then deter 
mined by extremizing the strain energy integral with respect to the family of 
admissable displacement states; the physically realizable state corresponds to 
the energy minimum. Quite logically, these concepts were immediately recog- 
nized as belonging to the analytical theory of self-adjoint systems, the vari- 
ational calculus and the theory of linear partial differential equations, as 
well as the corresponding principles in Hamiltonian mechanics inducing the 
Principle of Least Action and the Euler-Lagrange equations (cf. , Goldstein 
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ref 2) . Hence* the finite element method rapidly gained theoretical stat- 
ure as well as practical usefulness, albeit at the cost of conceptual restric- 
tion to linear problems in mechanics. 

In the mi d-1 960's, members of the engineering mechanics technical commu- 
nity, especially those with a background in mathematics or continuum mechanics 
but generally outside of structural mechanics, began to take note of the dem- 
onstrated power and versatility of the finite element procedure for analyzing 
complex problems. Particularly impressive was the ability of the algorithm 
to readily impose non- trivial (gradient) boundary condition constraints on 
irregularly shaped geometries, and to employ non-uniform discretizations of 
the problem solution domain (cf. Zienkiewicz, ref. 3). Each of these areas 
posed particular problems for more conventional numerical analysis procedures. 
Their apparent total alleviation by finite elements immediately prompted its 
extension to analyses of other linear field problems (cf. Zienkiewicz and 
Cheung, ref. 4) including steady heat conduction, and importantly subsonic 
potential flow. In this latter area, of our particular interest as a starting 
point, the differential equation is simply the linear Laplacian on the velocity 
perturbation potential function. However, the boundary conditions become 
specified on the behavior of its derivative normal to any surface, e.g., an 
airfoil, which in general is not parallel to a coordinate surface of the dif- 
ferential equation description. Because of the linearity, an energy-functional 
equivalent of the differential equation is readily established, including 
explicitly the non-trivial boundary condition statement. Using direct adapta- 
tion of structural computer code concepts, deVries and Norrie (ref. 5) obtained 
representative finite element solutions to two-dimensional aerodynamic flow in 
cascades using straightsided triangular elements. Computations for other two- 
dimensional aerodynamic configurations are given by Meissner (ref. 6); 

Sarpkaya and Hiriart (ref. 7) present results for an axisymmetric geometry 
involving a free surface of unknown location. Results using higher order 
finite element functionals, and curved-sided elements for improved geometric 
simulation, for two-dimensional potential flows are given by Thompson (ref. 8), 
Isaacs (ref. 9), and Hirsch and Warzee (ref. 10). The case for a specified 
freestream vorticity, "frozen" into an aerodynamic flow, is discussed by 
Vooren and Laburjere (ref. 11) using linear two-dimensional finite elements. 
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In all cases, the cited analyses correspond to the inviscid, steady, 
isoenergetic two-dimensional flow of an irrotational (or circulation preserv- 
ing), incompressible fluid, i ,e., potential flow aerodynamics. However, current 
high-lift aerodynamic geometries employ complex combinations of flap and slat 
systems. Powered lift configurations, using engine exhaust and/or mass ejector 
systems, are areas of currant aerodynamic research and development. All these 
flow fields depart significantly from linear potential flow, and viscous and/or 
turbulence effects form a large and important influence. Fortunately, interest 
in applying the finite element method of analysis to specifically non-linear 
flow fields was spawned in the late 1960's as well, with early attention 
focused on establishing an alternative theoretical foundation (since these 
problems did not appear self-adjoint!) The Method of Weighted Residuals was 
rediscovered (see Finlayson and Scriven, ref. 12), and emerged as a theoretical 
foundation for deriving finite element solution algorithms for arbitrarily non- 
linear partial differential equation systems. The linear potential flow case 
then became a special subclass of the more general formulation. Application 
to non-steady, two-dimensional aerodynamics is discussed by Bratanow and Ecer 
(ref. 13). Baker (ref. 14-16) presents numerous applications to a wide range 
of problems involving two- and three-dimensional viscous flows with turbulence 
and chemical reaction. Chan et al. (ref. 17) document extension of the 
finite element method to predictions in transonic aerodynamics. 

These computational results provide demonstration (by parts) of the poten- 
tial capabilities of the finite element method for analysis of complex aero- 
dynamic systems. However, for the "computational wind tunnel" to become a 
truly viable alternative to extensive physical testing, it is imperative that 
the computer conduct a compleat analysis. To do so, the geometrical modeling 
must be performed with high fidelity, accounting of viscous/inviscid inter- 
action must be intrinsic, and the inherent non-linearities of fluid mechanics 
must be accurately and economically captured by the computational mathematics. 
The two most important outputs from a computer study in aerodynamics are sur- 
face pressure distributions and overall drag; each can be accurately determined 
only by an adequate accounting for the entire physical system. Therefore, it 
is important to quantify the various flow disciplines, and to proceed through 
evaluation of the important analytical subsystems that constitute characteri- 
zation of a realistic subsonic aerodynamic configuration. 
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This report documents results of a study conducted for this purpose on 
the subject NASA Grant under the sponsorship of the Low Speed Aerodynamics 
Branch, STAD, NASA Langley Research Center. By resorting to the fundamental 
mathematics, the differential equation systems appropriate for description of 
various aerodynamic flow regimes are derived. Their basic structure is noted 
the underlying mathematical uniformity is utilized to subsequently establish 
the finite element solution algorithm applicable to all equation systems. The 
COMOC finite element computer program system, under development for analysis 
of problems in fluid and continuum mechanics, is briefly described. Subse- 
quently, numerical results on finite element analysis of the various problem 
classes of impact in aerodynamics are presented, with discussion of factors 
affecting solution accuracy, speed and adequacy as well as user-orientation 
of the analysis procedures. A summary of results completes the report. 

The work was primarily conducted while the senior author was on leave 
from Bell Aerospace Division of Textron as Visiting Associate Professor of 
Engineering Mechanics at Old Dominion University. We wish to acknowledge the 
long-term support given to research in finite element methods in fluid and 
continuum mechanics by Bell Aerospace, and the opportunity presented by Old 
Dominion University to carry out this work. 


NOMENCLATURE 

a boundary condition coefficient; unit vector 
A area 

b coefficient; unit vector 

B body force 

c coefficient; speed of sound 

Cp specific heat 

C chord; coefficient 

Cf skin friction 

D binary diffusion coefficient 

e alternating tensor 

E symmetric velocity gradient 

Ec Eckert Number 
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f function of known argument 

Fr Froude Numbe*' 

\ 

g function of known argument 

h static enthalpy; convection heat transfer coefficient 
H stagnation enthalpy 

i,j,k unit vectors of rectangular Cartesian coordinate system 
k thermal conductivity; constant 

K generalized diffusion coefficient 

l differential operator; turbulent mixing length 

L characteristic length; differential operator 

Le Lewis Number 

M Mach number; number of finite elements 

j n unit normal vector; coordinate normal to a curve 

j p pressure 

i Pr Prandtl Number 

I 

f q generalized dependent variable 

i Q heat addition; generalized discretized dependent variable 

j R perfect gas constant; domain of elliptic operator 

Re Reynolds Number 

| s coordinate parallel to a curve 

| S finite element assembly operator; boundary layer energy parameter 

I Sc Schmidt Number 

j t time 

I T temperature 

| u velocity vector 

j 

l U reference velocity 

i v perturbation velocity vector*, normal velocity 

x-j Cartesian coordinate system 

i! x # y„z rectangular coordinate system 

j; Y species mass fraction 

|| $ pressure gradient parameter 

j y ratio of specific heats 

j i 

|| SR closure of elliptic solution domain 

ji 6 boundary layer thickness; Kronecker delta 

j \ 

i 5 
f : 
i . 

I ; 
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A increment 

e kinematic turbulent viscosity; amplitude* belongs to 
k coefficient 

X multiplier; turbulence sublayer constant; period 
p molecular viscosity 

p density 

a integral kernel 

u stress tensor; integral kernel 

perturbation potential function 
$ scalar potential function; finite element functional 
X domain of initial -value operator 

X3 scalar component of vector potential 
f vector potential function 

to X3 scalar component cf vorticity; Van Driest coefficient 
fl vorticity vector; global solution domain 

{} column matrix 

[] square matrix 

U union 

fl intersection 

Superscripts 

T matrix transpose e 

a species identification - i,j # k 9 £ 

~ unit vector a 

■* vector m 

* approximate solution 0 

constrained, to closure p 

-L component normal to curve t 

// component parallel to curve u 

w 


Subscripts 

local reference condition 
tensor indices 
evaluated on lower surface 
mth finite element subdomain 
initial condition 
pressure coefficient 
partial derivative by time 
evaluated on upper surface 
evaluated at the wall 
global reference condition 
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THE FIELD EQUATIONS FOR AERODYNAMIC FLOW 


The complete description of a state point in multi-dimensional aero- 
dynamics is contained within the solution of a system of coupled, nonlinear, 
second order partial differential equations describing the conservation of 
mass, linear (or angular) momentum, and energy. To establish unique solutions, 
closure of the system is required by specification of appropriate constitutive 
behavior and boundary conditions. In conventional vector notation, and then 
repeated in the preferred Cartesian tensor notation (where the subscript comma 
implies the gradient operator, subscript comma t is the time partial de- 
rivative, and the subscript semicolon implies the generalized divergence 
operator, see ref. 15), the conservation form of the differential equation 
system is 


r= 

3t 1 


[pul 


—■(pu) = -V* [pUU -■?] + pfr 


(1A) 


(2A) 


I pH - p] - -V-[pu H -t-u - kVT] + pq 


3t 


(pY a )= -V-IpuY 01 - pDVY a ] 


(3A) 

(4A) 


P, t = -tpy.1.. ■ 

( pu i ) *t = “ EpU i u j + P 6 ij " T ij 3 ;j + P 8 ! 
(pH - - T.jUj - kT 9 ^ ] . . +.pQ 

(pY a ) »t “ "tP u 1 Y a - pDY*,.].. 


( 1 ) 

( 2 ) 

(3) 

(4) 


The dependent variables in eqn.(l)-(4) have their usual interpretation in 
fluid mechanics. The mass flow vector is pu. where p is the mass density 
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and u.j is the local velocity vector. In eqn.(2), EL is an appropriate 
body force s and t.. is the diffusive (viscous) stress tensor. £qn.(3) is 

* J 

written on stagnation enthalpy, H,p is the static pressures T is the 
static temperature, and Q is the local heat generation rate. For multiple 
species flows, as might occur for example in powered lift configurations, Y & 
is the mass fraction of the species in eqn.(4), and D is the binary 

diffusion coefficient. Eqn.(4) also provides the means for tracing the trans- 
port of distinct flow field components. 

As stated, the solution of eqn.(l}-(4) requires specification of consti- 
tutive relationships between the dependent variables and the diffusion coef- 
ficient D , the stress tensor t.., and the effective viscosity and thermal 

* J 

conductivity. For compressible flows, an equation of state relating the 
thermodynamic variables is also required. Since predictions in both laminar 
and turbulent flows are required, and since the time-averaged form of eqn.(l)- 
(4) appears identical to those presented, by appropriate interpretation of the 
stress tensor, the general form of closure can be written as 

T^-tulD+peymu^j + Uj.,) (5) 


In eqn.(5), p identifies the laminar molecular viscosity which is tempera- 
ture dependent, and is the effective transport tensor due to turbulence. 

As indicated, the diffusive stress tensor is assumed a functional of the sym- 
metric velocity gradient (cf., Donaldson et al., ref. 18). For strictly lam- 
iner flow of a Newtonian fluid, eqn.(5) embodies Stokes' viscosity law 



= 2ij[E ij4 Viji 


( 6 ) 


which displays the linear functional dependence on the symmetric velocity 
gradient as 


_ 1 




( 7 ) 


The thermodynamic properties are typically expressed as 
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( 8 ) 


p - p(p»V a »T) 



(9) 


D a * D(Y a fl u i ) 


£ 10 ) 


Rather than attempting direct solution of eqn.(l)-(10) , it is convenient 
to non-dimensional ize all variables to extract the useful non-dimensional 
groupings. Using standard procedures (ref. 15.), the non-dimensional form of 
eqn.(l}"(4) is 



(11) 


H>,t = - [p°i u j + j " w T i j J SJ + 
(pH - (Ec)p), t = -[pu, H - § Ty u. - jfev H 

(p v “!’t = - h ya - W ^ V “ ’i] ;i 



_J 

^ Pr^, 


pQ 



( 12 ) 

(13) 

(14) 


In eqn. (li)-(14) , the important non-dimensional parameters for fluid mechanics 
are identified as 
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Reynolds Number: 


05 ) 


Prandtl Number: 
Eckert Number: 
Froude Number: 

Lewis Number: 


pJJ L 
Re a-21=- 

\i 

* m 


C D^ 

py» = . P. 
rr _ k 


Ec = 


U 2 

CD 

C p T co 

^ CO 

u 2 


Fr = 


Lg 


Le = pp 

U 


( 15 ) 

(17) 

08 ) 

09 ) 


It is also useful to define the Schmidt number from the above as 

C - Pr 
Sc = Le Re 

Note that, for a thermodynamically perfect fluid 

Ec = (y-DM^ (20) 

where y is the ratio of specific heats, and is the reference Mach 
number defined as 

U 

M = — — 

w (2D 

where R is the universal gas constant divided by the molecular weight. 

The solution of ecjn.(n)-(14) is a formidable task, and there are several 
alternative forms and simplifications to the system that can render such solu- 
tion more straightforward for non-trivial problem classes in aerodynamics. 

The Navier-Stokes equations primarily describe the behavior of the mass flow 
vector field, pu.. From vector field theory, we recognize that without loss 
of generality, we can describe pu^ in terms of a scalar and a vector poten- 


12 



tial function, <f> and $ respectively, as 

pij = -p I + ^ x f 


(22A) 


pu- - -p$,^ + e. 


ijk ^ksj 


( 22 ) 


The first equation employs conventional vector notation; the second illus- 
trates the preferred indical form, and the decomposition in eqn.(22) is unique 
to within an arbitrary constant, provided $ satisfies the "gauge" condition, 
i.e., W = “ 0 (cf. , 19, App. I.). In terms common to fluid mech- 

anicSj $ is called the (total) potential function, and $ is the three- 
dimensional equivalent of the familiar two-dimensional stream function. 

Viewing eqn.(ll), we observe that for steady or incompressible flows, the 
massflow vector field also satisfies the "gauge" condition, i.e., it is diver- 
gence-free, Inserting eqn.(22) into (11), and from the skew-symmetry proper- 
ties of the Cartesian alternator, we obtain 

( pU i ) ;i = Cp$s i J ;i = 0 (23) 

An elementary solution to eqn.(23) is that $ vanishes everywhere. Hence, 
for all cases where eqn.(ll) describes steady or incompressible flow, we may 
choose to analyze the flow using the transformation 

pu i = e ijk f k,j (24) 

In this case, pu^ is said to be a solenoidal vector field. An important 
additional aspect of this problem class is that the flow may be explicitly 
rotational. The vorticity fis is defined as the curl of velocity field; hence 

fi i 5 e tilc u k:i (25) 
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Substituting eqn. (24) into (25), utilizing the symmetry properties of alter- 
nator contractions, and noting that ^ satisfies the gauge condition, the 
compatibility relation between stream function and vorticity is directly ob- 
tained as 


£ 2 , 



(26) 


Eqn. (26) may find considerable use in analysis of certain problem classes in 
aerodynami cs . 

We have an alternative choice for the aerodynami ca 11 y important case of 
irrotational flow, wherein the curl of eqn. (22), or eqn. (25) vanishes identi- 
cally by assumption. Again using the skew-symmetric properties of the alter- 
nator, and for satisfying the gauge condition, one obtains from eqn. (22) 


e jik ( pu i ^ k - 0 ^k;jj 


(27) 


If we set l l' equal to zero everywhere, for irrotational flows we then have 
the identity. 


PU| = -p§,* 


(28) 


In this case, pu^ is said to be a lamellar vector field. 

We may elect to employ the primitive variable description for pu- , or 
may selectively employ the definitions in eqn. (24), (25) and/or (28). The 
particular choice depends upon the features of the differential equation sys- 
tem derived from eqr.. ( 1 1 ) - ( 1 4) for each case of aerodynamic value. 


Aerodynamic Potential Flow 

The steady, isentropic flow of an ir.viscid, single-species perfect gas 
has been the focal point of research in aerodynamics for well over 70 years. 

A specific accounting of the simplifying assumptions is readily accomplished 
for the parent differential equation system, eqn.(ll)-(14). Neglecting body 
forces, and discarding the species continuity equation as redundant, for this 
problem class we have 
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V, 


(29) 


(pu n -) ;1 = 0 


H u j>;i * -P»j 


(30) 


c 2 p #i = p #1 


(31) 


Since isentropic implies adiabatic and reversible, eqn. (13) is identically 
satisfied by H s constant . Eqn.(31), an alternative expression for this 

form of eqn.{13), defines the speed of sound in an isoenergetic perfect gas. 
Since the flow is assumed inviscid, the dispersive stress tensor, eqn. (5), 
vanishes identically. Eqn.(29)-(31) can be conveniently combined into a 
single differential equation on the velocity field, u^ . Substituting 
eqn.(31) plus the expansion of eqn.(29) into eqn.(3Q), one directly obtains 
the differential constraint for steady, inviscid isentropic flow as 





= o 


(32) 


In eqn. (32) , recall that c is the local speed of sound, conveniently defined 
in terms of a stagnation reference - condi ti on and the local velocity field by 
the isentropic energy equation in the form 


c 


2 




(33) 


Eqn. (32) is a highly non-linear, first order partial differential equa- 
tion written on the field behavior of the velocity vector u^ . It is valid 
for all Mach numbers, and as a function of local Mach number may selectively 
display elliptic, parabolic, and hyperbolic differential character. An impor- 
tant subclass corresponds to the additional constraint that the flow be irro- 
tational. For this case, a useful restatement of eqn. (32) is obtained using 
eqn. (28). By direct substitution, the second order partial differential equa- 
tion for determination of the distribution of the scalar potential function $ 
is 


15 



Eqn . (34) displays the mixed differential character as well. For i = j , the 

term in the bracket becomes of the form [1 - » where M^.j is the 

local Mach number of the flow in the coordinate direction Hence, for 

subsonic flows, eqn. (34) is elliptic, and from the theory of the solution of 

partial differential equations we require knowledge of an algebraic combina- 
tion of $ and its normal derivative, 5,^ around the complete closure of 
the solution domain. For supersonic flows, one of the bracketed terms may 
become negative, and eqn. (34) then displays a mixed elliptic-hyperbolic char- 
acter. In this instance, boundary specifications on <p and/or its normal 
derivative are required on some non-characteristic curve. In the intermediate 
(transonic) speed range, eqn. (34) can exhibit a globally elliptic character 
within embedded hyperbolic and/or parabolic sub-domains, and the numerical 
solution procedure must recognize these distinct regions. 

Regarding boundary conditions, we observe that specification of velocity 
on the closure of the solution domain will correspond to a gradient boundary 
condition specification on the scalar potential function, see eqn. (28). On an 
arbitrarily oriented closure segment, identified by a unit outward pointing 
normal vector, n^ , eqn. (28) states that 

H \ - -Vk 8 - ul < 35) 

.1 

where J is the scalar component of velocity parallel to . As a special 
case, eqn. (35) states that the normal derivative of $ must vanish on an 
impervious aerodynamic surface, sinoe thereupon the local velocity vector must 
lie tangent to the surface. 

The numerical solution of the general non-linear form of eqn. (34) is not 
commonly attempted. For specification of additional simplifying assumptions, 
which constrains the generality of the flow field description, alternative 
forms for eqn. (34) that are more tractable for numerical solutions using con- 
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ventional techniques can be obtained. By and large, the fundamental step is 
re-definition of the scalar potential function, eqn.(28), under the assumption 
of the preponderant existence of a preferred flow direction. Assuming this 
direction aligned with the x-j coordinate axis, and identifying a reference 
velocity (Uj as parallel to this direction, re-define the local velocity 
field as 


pu i = pU w 6 in + pv. (36) 

Under the assumption that the perturbation velocity component, v* , is small 
in magnitude in comparison to the freestream reference velocity, we can define 
a corresponding perturbation potential function as 

pv i = -p<f> M (37) 

Then, eqn.(28) takes the form 

pu* = -p$ #1 . = pU TO 6 • 1 -p<f) si (38) 

Substituting eqn.(38) into eqn.(32), and altering the reference condition in 
eqn.(33) to that corresponding to U ro , and proceeding through the well known 
order of magnitude analysis (cf. , ref. 20), for the subsonic and supersonic 
slender body approximation one obtains the linearized form of eqn.(34) 
written on either the total or perturbation potential function as 

1 - + *;22 + ®;33 = 0 (39) 

The boundary conditions for solution of eqn.(39) are obtained from the 
definition in eqn.(38) in the manner analagous to the specification in 
eqn.\35). Eqn.(39) displays the mixed elliptic-parabolic-hyperbolic character 

of the full non-linear form, eqn.(34). For the final additional assumption of 
zero Mach number, or equivalently incompressible flow, eqn.(39) becomes the 
Laplacian on § , which is a linear elliptic-boundary value problem. It is 
this form whose numerical solution is typically computed for subsonic flows. 


17 


in combination with an independent variable transformation on the coordi- 
nate axis to account for a finite subsonic Mach number. 

Flow fields of the subject problem class are also amenable to analysis 
using the definition of stream function, eqn.{24), provided the parent flow is 
either incompressible or steady. For the additional assumption of irrotational, 
eqn.(26) must vanish identically from the definition in eqn.(25); hence 

[p*k;l]., "(?),, ' t ’i;k = 0 (4°) 

Eqn.(40) states that the vector potential function, ^ , satisfies a second- 
order linear partial differential ecuation of the Laplacian type provided 
that the stream function satisfies the gauge condition. Since, from its def- 
inition, density can never become negative, and since the first-order differ- 
ential term in eqn.(40) does not affect overall characterization, eqn,(40) is 
uniformly elliptic and its solution requires knowledge of and/or the 
normal derivative, T..n. on the complete closure of the solution domain. 

K s J J 

The vector potential function is simply a generalized concept of the two- 
dimensional stream function familiar to all. Therefore, equal to a con- 
stant implies existence of a stream hypersurface across which mass flow van- 
ishes. Forming the outer product of the defining equation, eqn.(24), with an 
arbitrary unit vector yields 

P V* = *Hjk '•'kj < 41 > 

For a £ parallel to the velocity vector u. , the specified mass flux across a 
segment of the solution domain closure determines the distribution of ^ on 
the corresponding surface. Conversely, for lying parallel to the velocity 
vector u.j , eqn.(41) yields a normal gradient-type boundary condition con- 
straint on stream function in terms of the parallel velocity component of the 
form 


.n. = -U," 
P J k 


(42) 
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It is instructive to contrast the form of the boundary condition constraints 
for the two potential functions, see eqn.{35) and (42). As a final simplifi- 
cation, in the instance of constant density flow, note that eqn.(40) degene- 
rates to the linear Laplacian written on ^ . In all instances, since both 
the scalar and vector potential functions are defined only to within an arbi- 
trary constant, see eqn.(22), each function may be appropriately specified at 
some point on the solution domain closure for convenience. 

Analysis of rotational, inviscid isentropic flow fields can be accomp- 
lished using either the scalar or vector potential function as well. Regard- 
ing the former, for linearized potential flow theory, the Kutta-Joukowski 
hypothesis states that the correct rotational flow pattern of any family of 
flows is the one flow with a finite velocity at the trailing edge of an air- 
foil. Hence, from the concept governing two-dimensional, irrotational incom- 
pressible flow about bodies resembling subsonic airfoils, one can complete an 
analysis of a family of flows for a given boundary and a given freestream 
velocity. These flows differ primarily in the amount of circulation; however, 
all but one will have an infinite velocity at the trailing edge. With use of 
Kutta condition, use of eqn.(34) or any of its simplifications is not contra- 
indicated for rotational inviscid flows. 

For rotational incompressible or steady flows, use of the vector poten- 
tial function is straightforward. The derived compatibility eqn.(26) becomes 
in this instance 




P k,j 


= -SI. 




(43) 


Eqn.(43) states that the distribution of the stream function will be altered 
by the imbedded rotational character of the flow. Provided th$ flow is steady 
and inviscid, an analysis (discussed in the next section) shows that the 
rotational freestream will be "frozen" into the computed streamline distribu- 
tion about an aerodynamic shape. Hence, by establishing an iterative pro- 
cedure for solution of eqn.(43), one theoretically can determine the rota- 
tional flow about an arbitrarily shaped aerodynamic surface, with embedded 
freestream vorticity, using the vector potential function. 
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Aerodynamic Viscous Flow 

The complete analytical characterization of an aerodynamic system requires 
an accounting for viscous effects. For conventional aerodynamic shapes at 
small angle of attack, these effects are essentially confined to the thin layer 
of fluid in direct contact with the aerodynamic surface. The remainder of the 
flow field is essentially free of viscosity effects and completely amenable to 
analysis using inviscid flow field concepts. These thin viscous layers are 
called boundary layer flows, and their analysis and understanding constitutes 
an extremely important branch in engineering and aerodynamics. Under a first- 
order of magnitude analysis (cf. , Schlichting, ref. 21), for the steady bound- 
ary layer flow of a viscous compressible fluid described by the velocity 
vector 

u. - u 1 e ] + u 2 c 2 * u 3 e 3 (44) 

where the unit vector triad e.. lies parallel to coordinate curves of a 
general Cartesian curvilinear coordinate system, the three-dimensional bound- 
ary layer equivalent of eqn.(ll) -(14) takes the form 


o = (pu^.. 

(45) 

pu i u l;i = Re [ ty * pe 12 )u l ;2]. 2 " p »l 

(46) 

0 « - P , 2 

(47) 

pu i u 3;i = Re [^ u + pe l 2^ u 3 52 ] . 2 “ p =*3 

(48) 

pu i H, i RePr + p£ 12) Hs 2 * Re ( u l;2 ' u 3;2) 

J ;2 

(49) 

pu i Y ’i = W * pe ig) ^ 2 ] . 2 

(50) 
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In the form presented, the use of a body-oriented orthogonal coordinate system 
attached to the aerodynamic surface is assumed, where x-j is parallel to the 
predominate direction of flow and x 2 is assumed to lie everywhere perpendic- 
ular to the aerodynamic surface. For aerodynamic surfaces with local curva- 
ture distributions, the generalized semicolon subscript (;} differential nota- 
tion will yield additional terms stemming from use of a curvilinear coordinate 
system in equations (45)-{50). Note that the flow may be laminar, transitional 
or turbulent. 

Perhaps the most significant feature of the boundary layer equations, and 
their solution, is contained in eqn.(47). Under the first-order of nagnitude 
simplifications, the transverse momentum equation yields the significant fact 
that the pressure distribution is constant i-?roughout the bounday layer thick- 
ness. Hence, in the remaining momentum eqn.(46) and (48), pressure appears as 
a parameter only, being obtained from the inviscid flow field solution. There- 
fore, eqn.(46) and (48) are solved respectively as initial value problems for 
u-j and u^ with the pressure replaced by p e , while eqn.(45) is recast into 
an initial -value problem for the x 2 distribution of the velocity component 
normal to the surface, u 2 . For non-isoenergetic flows, eqn.(50) is also 
solved as an initial value problem for the distribution of stagnation enthalpy, 
H . It should be noted that the conventional two-dimensional boundary layer 
equations are a sub-set of the presented form, obtained by deleting aqn.(48) 
and excluding 3 as an admissible index for summation in eqn.(45)-(5Q). 
Finally, for single species flow, eqn.(50) may be deleted completely. 

The initial and boundary conditions appropriate for solution of eqn.(45)- 
(50) are obtained in a straightforward manner, as all except (45) represent a 
general initial -valued, two-point elliptic boundary value problem in mathe- 
matics. On the aerodynamic surface, the no-slip boundary condition for the 
tangential velocity field yields 

u 1 (x-|,0,x 3 ) = u 3 (XpO,x 3 ) - 0 (51) 

To account for suction or blowing at the aerodynamic surface, the correspond- 
ing constraint for solution of eqn.(45) becomes 

U 2 U 1 — v(x 1 ,x 3 ) (52) 
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where v is the specified distribution of normal velocity. Similarly, at the 
frees tream juncture (x 2 = 6) between the boundary layer and potential flow, 
we require the boundary layer velocity distribution to be equal to that of the 


external inviscid flow, denoted by a subscript e , as 

u^Xpfi^) = u (x 1# x 3 ) (53) i 

Typically, solution of eqn.(50) at the surface involves a convection-type 
boundary condition written in terms of the local static temperature as 

-kT, i n^ = h[T - T r ] (54) 

At the freestream, the normal gradient of H typically vanishes as 

H(x.j ,S,x 3 ) , • n. = 0 (55) 

Finally, the form of eqn.(55) is typically appropriate at x 2 = 0 and x 2 - <5 
for solution of eqn.{51), with H replaced by Y a . 

Since eqn.(46), (48)- (51 ) are also initial value problems in the x-j 


coordinate direction, an initial profile distribution for each appropriate 
dependent variable is required. Denoting q as a generalized dependent 
variable, the solution of the equation system is initialized by specification 
of the form 


q(0 a x 2 ,x 3 ) — (**6) I 

j 

The presented two- and three-dimensional boundary layer equations are the j 

most tractible form of the parent Navier-Stokes equations for numerical solu- j 

tion. Consequently, their range of applicability is considerably limited by j 

geometric or flow field considerations. For merging viscous flows near trail- } 

i 

ing edges, or flap and slat combinations, or for powered lift configurations 1 

with their associated thick viscous non-isoenergetic flow fields, a more com- ! 

I 

prehensive form of differential equation description is required. For steady j 

three-dimensional flows of this type, in bounded or open domains and meeting 
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certain requirements, a significant simplification can be made to the parent 
Navier-Stokes equations that renders three-dimensional solutions considerably 
more retractable with present computer hardware. This approximation, now 
known as the “parabolic Navier-Stokes equations," describes steady three- 
dimensional flows wherein: 1) a predominant flow direction is uniformly 
present, 2} in this direction (only) diffusion processes are negligible com- 
pared with convection, and 3) no disturbances are propogated upstream anti- 
parallel to this direction. For the velocity vector identified in eqn.(44), 
and for the same coordinate description of the solution space, the three- 
dimensional parabolic Navier-Stokes equations are of the form 


0= (pu.).. 

O' 4 *) r 


pu j u i;j " Re 


* pe ijc) ^i ;k 


;k 


- p ; 


(57) 


(58) 


0 - *«) 
pu i H, i Ri 


(P + P^) 


L 


Pr 




;k 


1 - Pr 


( u + p e ik) 


Pr 


<u Mk 


;k 


Sc - Pr 


ScPr 


(<* + p-ik) E 


a 


;k 


(59) 


p Ul Y“. 


0 ■ \i) 


Re 


+ pe iO 

' »L 


Sc 


;k 


(60) 


The dominant differences between the parabolic Navier-Stokes equation 
system, (57) -(60), and the three-dimensional boundary layer eqn.(45)-(50), 
relate to dimensionality of the boundary value character and to the appearance 
of pressure. For boundary layer flow, eqn.(47) yielded a uniform pressure 
distribution imposed through the thickness of the boundary layer. In the 
parabolic Navier-Stokes system, pressure appears as a dependent variable and 


its solution is required. The several forms of the pressure algorithm, to be 
discussed, depend upon the particular geometrical configuration under study. 
The other significant feature of the parabolic Navier-Stokes system is that 
the diffusion term {first right side term modified by (l-S^J/Re) is two- 
dimensional, spanning the plane whose normal is parallel to the direction of 
predominant flow. (The subscript bar notation indicates no summation; for 
k-1 only, O-^i) vanishes.) Hence, boundary condition statements on the 
three components of velocity as well as enthalpy and mass fraction must be 
specified everywhere on the closure of this two-dimensional plane. Identify- 
ing q for a generalized dependent variable, most physically realistic bound- 
ary condition constraints are of the form 


a Cl) (x i )q + a^Cx^q,^ = a (3) U.) (61) 

In eqn.(61), the a^ are user-specified coefficients; note that the corres- 
ponding boundary layer equation statements, eqn. (51 )-(55) , are all special 
cases of eqn. (61). As before, since the parabolic Navier-Stokes equation 
system contains an initial value description for each dependent variable, 
specification of the form of eqn. (56) is required for each dependent variable. 

Dependent upon the geometry of the flow configuration, eqn. (57) and the 
three equations in (58) are selectively altered to obtain solution of the 
three components of velocity and the pressure distribution. For boundary 
layer-type flows (termed boundary region flow), since eqn. (57)-(60) encompass 
eqn.(45)-(50) as a special subset, the assumption of vanishing normal distri- 
bution of pressure may be valid. Hence, the inviscid pressure distribution vs 
imposed upon the flow field. Pressure is again decoupled from the solution, 
the second of eqn. (58) is discarded, and eqn. (57) is employed to solve for the 
x ? distribution of the corresponding velocity component u 2 . For solution 
domains totally bounded by solid walls, as might occur in inlets or exhaust 
ducts for example, the theory of Patankar and Spalding (ref. 22) can be used 
to achieve a pressure solution. Their theoretical procedure involves a split- 
ting of the pressure field computation by decomposition of the local static 
pressure into the form 
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p(x i ) = p(x-|) + p’(x 2 »x 3 } 


( 62 ) 


Determination of the x-j component, p , is obtained by integrating eqn.(57) 
oyer the cross-sectional area of the duct and accounting for the influences of 
wall shear (t w ), area change (A(Xq)), and heat or mass addition to yield the 
solution form 

2^* f(pu r A(x,J. V A t) (63) 

Note that eqn.(63) is an ordinary differential equation; its solution deter- 
mines the (assumed uniform) axial pressure gradient appropriate for solution 
of the first of eqn.(58). To obtain solution for the distribution, in the 
transverse plane p'fxgsX^) , the divergence of the second two equations of 
(58) is taken and coupled with eqn.(57) to yield the solution form 


O - 6y)p'. kk * f (pu r V“, T) (64) 

As written, eqn.(64) is a two-dimensional elliptic boundary value problem of 
the Poisson type; the right side is a specified function of its arguments. 

The previous approach is inappropriate for flows in completely unbounded 
domains. In this instance and coupled with an equation of state, aqn. ( 57) can 
be cast as a pure initial value problem on the three-dimensional pressure 
distribution, p(x.) , In this case, eqn.(57) takes the form 


dp(x.) 

dx-j 




+ n - « kl ) (pu k ;. k 


(65) 


where we have assumed valid the perfect gas law. Note that eqn.(65) is 
totally inappropriate for confined flows, since the right side becomes infin- 
ite at a wall where Uq vanishes. For three-dimensional flows bounded by an 
aerodynamic surface and an inviscid freestream, the second of eqn..(58) yields 
an appropriate specific form for determination of transverse pressure distri- 
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bution as 


dp(x 2 ) 

dX n 


= f(pu. 


T, P e ) 


( 66 ) 


For the alternative case of a viscous flow field imbedded between two poten- 
tial flows, as might occur in the trailing edge region downstream of an air- 
foil for example, a recasting of eqn.(57) into the form of eqn. (64) using the 
concepts of eqr.(62) may be required. In all instances, the computation of 
pressure typically involves application of initial- value techniques coupled 
with the explicit assumption that pressure variation is a local phenomenon 
unaffected by downstream influences. For flow fields where this assumption is 
violated, use of the parabolic Navier-Stokes equation systems is probably 
contraindicated. 

Solution of the boundary layer or parabolic forms of the Navier-Stokes 
equations is appropriate for flow fields where curvature effects are suffi- 
ciently modest such that streamwise separation does not occur. However, at 
larger angles of attack for an airfoil for example, the created adverse 
pressure gradient will retard the parent uni di recti onal-type flow to the point 
where the streamwise momentum is insufficient to keep the flow attached to the 
airfoil. The phenomenon of separation occurs, whereby the remaining unidirec- 
tional aerodynamic flow is separated from the aerodynamic surface by a suit- 
ably sized region of highly rotational viscous flow. For fully three-dimen- 
sional geometric configurations, characterization of these flows requires 
solution of the full Navier-Stokes eqn. {11 )-(14). However, since solution of 
this form is typically not tractable on current generation computers, full 
Navier-Stokes solutions are generally limited to two-dimensional configurations. 

The form presented as eqn. (11)- (14) can be used for solution of two- 
dimensional viscous flows. However, for steady or incompressible cases, an 
alternative formulation can be completed that takes advantage of the vector 
field character of eqn. (11). For this case, eqn. (24) takes the specific form 


pU i = e 3ij *3,j 


= e 3ij 


(67) 
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It is significant to note that only the Xg scalar- component of the vector 
potential function is required to characterize two-dimensional flows. Similar- 
ly, eqn.(25) becomes 


fi 3 = to = e 3ij Uj ^ (68) 

The compatibility equation, eqn. (26), becomes for two-dimensional flow, using 
eqn.(67) and (68) 



In terms of stream function and vorticity, the momentum equation (12), of 
the Navier-Stokes system for two-dimensional steady or incompressible flows, 
becomes considerably simplified (ref. 15) to the form of the vorticity trans- 
port equation. 



The most significant feature of eqn.(70) is the disappearance of pressure as a 
coupled dynamic variable. Hence, eqn.(69) and (70) constitute a closed system 
for determination of two-dimensional rotational viscous aerodynamic flows. 

The stagnation enthalpy and species continuity equations, eqn.(13)-(14) ,, 
remain as presented except for substitution of eqn.(67) for the terms involv- 
ing velocity. 

As stated before, eqn.(69) is an elliptic boundary value problem since 
density never vanishes and is always positive. Therefore, for boundary con- 
ditions we require knowledge of \jj and/or its normal derivative everywhere 
on the closure of the solution domain. For finite Reynolds number, eqn.(70) 
is also an elliptic boundary value problem coupled with initial -value behavior 
due to the time derivative (for incompressible flows only). Hence, we also 
require knowledge of w or its normal derivative everywhere on the closure 
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of the solution domain. Since the defining eqn.(68) is valid throughout the 
solution domain as well as its closure, boundary condition specifications on 
vorticity are obtained from eqn . (68) and (69). The numerically consequential 
vorticity boundary condition becomes the imposition of the no-slip wall; the 
equivalent embodiment in vorticity is 

(71) 

an 2 

In eqn. (71), n is the coordinate normal to the local closure segment of the 
solution domain. To evaluate eqn. (71), it is necessary to form the second 
derivative of the streamfunction distribution on the closure. Several forms 
have been determined as appropriate, see Roache (ref. 23, Section III.C). 

The pressure distribution may be recovered from the vorti city-stream 
function characterization of a Navier-Stokes solution by evaluating the momen- 
tum equation, eqn. (12). Two alternatives exist; eqn. (12) can be differentiated 
by the divergence operator yielding a boundary-value specification on pressure 
involving the Laplacian operator. However, since we are primarily interested 
in pressure distributions on aerodynamic surfaces, an alternative formulation 
exists which can be of value. Form the vector contraction of eqn. (12) with an 
infinitisimal displacement vector dx.. and integrate over any curve; this 
yields 

/ P*4 a 1J dx 1 = ■/ [ p Vj - IS T ij ] pu i d *i (72) 

f Mote that the left side of eqn. (72) is the integral of a perfect differential 

and thus independent of path. Hence, the pressure at any point in the field 
can be determined, in comparison to some reference value, by integration over 
j an arbitrary (the easiest) path between the two points. Denoting the inte- 

\ grals of perfect differentials by A , completing the indicated integrals, and 

!; performing numerous integration by parts, the final form for eqn. (72) becomes 

I. 

I. - ' 
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In eqn.(73), the first curly bracket contains terms which are point dependent, 
the result of integrating perfect differentials. The second and third brack- 
ets require integration, thus selection of path. While the appearance of 
eqn.(73) is formidable, its numerical evaluation utilizes well-known techniques. 

As a final observation, we noted in our earlier discussions that for 
steady, constant- density inviscid rotational flows, the vorticity was frozen 
into the streamline distribution. This can be readily determined from eqn.(70), 
the vorticity transport equation. The complete first term vanishes identically 
for inviscid flow, while the left and last right side terms vanish identically 
for steady, constant density flows. Hence, the sole term remaining for this 
case is due to convection. In terms of a streamline coordinate, denoted as 
s(x.), eqn.(70) can be written in scalar notation in terms of the velocity as 



(74) 


Since the streamline speed u does not vanish, satisfaction of eqn.(74) is 
obtained only for vorticity cu vanishing identically or not varying in the 
streamline coordinate direction. Hence, the previous observation for aero- 
dynamic potential flow is confirmed from the vorticity transport equation. It 
is useful to note that the form of eqn.(74) is quite identical to Bernoulli's 
law for determination of pressure distributions along inviscid streamlines. 


uf.l?p (75) 

as p H 

Eqn . (75) is utilized for determination of surface pressure distributions using 
the velocity distribution computed from an inviscid flow field analysis. 
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FINITE ELEMENT SOLUTION ALGORITHM 


Various forms of the Navier-Stokes equations have been developed as 
appropriate for analysis of distinct configurations in aerodynamic flow. For 
an inviscid irrotational analysis, eqn.(34) is the preferred parent form, 
written on scalar potential function, with eqn.(39) obtained following some 
linearizing assumptions. Vector streamfunction is applicable to rotational 
flows as well; its distribution is established via solution of eqn.(40) or (43) 
appropriately. Viscous flow analyses have fallen into several categories. For 
a completely attached viscous boundary layer between the aerodynamic surface 
and the freestream, solution of the two- or three-dimensional boundary layer 
eqn.(45)-(50) is appropriate. For merging viscous flows near trailing edges, 
or flap and slat combinations, or for powered lift configurations with thick 
viscous non-isoenergetic flow fields, the more comprehensive parabolic Navier- 
Stokes system is appropriate, eqn. (57) -(60) , provided the flow remains unidi- 
rectional. Finally, analysis of omnidirectional viscous flows requires solu- 
tion of the complete Navier-Stokes equations, either in the parent form of 
eqn. {11)-(14) or for two-dimensional configurations, eqn.(69)-(70). For the 
subsonic flows of present interest, each member of the derived differential 
equation systems (except for instances with specific retention of the contin- 
uity eqn. (11)) is uniformly cast as an elliptic boundary value problem of math- 
ematical physics with individual instances of initial -value coupling. Specif- 
ically, identifying q(x.j,x) as a generalized dependent variable of interest, 
each of the subject partial differential equations is a special case of the 
general, second-order non-linear partial differential equation 


L(q) = k! 




+ ftq.q^x.) - g(g,x) s 0 




(76) 


In eqn. (76), f and g are specified functions of their arguments, x is 
identified with x-j for boundary layer or parabolic flows or time for trans- 
ient flows, and the x.. are the coordinates for which' second- order derivatives 
exist in the lead term. For this term, k is a scalar constant and k..(q) 

* J 

is the generalized diffusion tensor, and both become uniquely specified by 
identification of q with each of the dependent variables associated with the 
particular equation system of interest. The finite element solution algorithm 
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is based upon the assumption that L(q) is uniformly parabolic within a 
bounded open domain £2 ; that is, the lead term in eqn.(76) is uniformly 
elliptic within its domain R , with closure 3R , where 

Si = R * (Xq»x) (77) 

and x 0 < X . If eqn. (76) is uniformly parabolic, unique solutions for q 
are obtained upon specification of functional constraints on the closure of SI, 
3ft = 3R x[x 0 »x) » and an initial -condition specification on RU3R * x 0 • For 
constraints on 3S2 , the general form relates the function and its normal 
derivative everywhere on the closure 3R as 

&(q) = a^ q(x!. »x) + qCx.. ,x) ^ n^ ~ a ^ = 0 (78) 

In eqn.(78), the {x,.,x} are user-specified coefficients, the superscript 

bar notation constrains x- to 3R , and n. is the local outward-pointing 

■ J 

unit normal vector. For an initial distribution, assume that 

q(x. ,x 0 ) = %( x i) (79) 

is given throughout RU3R x X 0 * 

The finite element solution algorithm is established for the equation 
system (76)-(79) by using the method of weighted residuals (MWR) formulated on 
a local basis. Since eqn. (76) is ‘valid throughout n , it is valid within 
disjoint interior subdomains described by (x^ sX) eR m x [Xq»x) * called 
finite elements, wherein UR^ = R . An approximate solution for q within 
Rjn xtXgsX) > called qj( x ^x) s is formed by expansion into a series solution 
of the form 

q * (Xi,x) ={«(x i )} T {q(x)} m (80) 

In eqn. (80), the functionals ^(x^) are subsets of a function set that is 

on coefficients Q^(x) represent the unknown 
X-dependent values of q*{x^,x) at specific locations interior to R m and 
on the closure 3R Jn , called nodes of the finite element discretization of R . 
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complete on R^. The expans i 


To establish the values taken by the expansion coefficients in eqn.(80), 
require that the local error in the approximate solution to both the differ- 
ential equation L(q*) and the boundary condition statement £(q*) n for 
aiyiBR f 0 , be rendered orthogonal to the space of the approximation func- 
tions. By employing an algebraic multiplier A, the resultant equation sets 
can be combined as 

= {0} (81) 

where jj is the mapping function from the finite element subspace R m to the 
global domain R , commonly termed the assembly operator. The number of 
eqn.(81) prior to assembly is identical with the number of node points of the 
finite element R [n . 

Eqn„(81) forms the basic operation of the finite element solution algo- 
rithm and of the COMOC computer program to be described. The lead term can be 
rearranged, and A determined by means of a Green-Gauss theorem: 

{*( x i>} K [ K ij q m’i].. dT = K f, R K ij q m’i"j da 

S J 

' K I ( t(x i ) } ‘i K ij q mM dT (82) 

K m 



/ (Hx,)} L(q*)dT-xJ A (q*) 

L \ V 3R 


da 


For 3Rn3R m nonvanishing in eqn.(82), the corresponding segment of the closed- 
surface integral will cancel the boundary condition contribution (eqn.(81)) 
by identifying Xa^ with k of eqn.(76). The contributions to the closed- 
surface integral eqn.(82), where 3R m n3R = 0 , can be made to vanish (ref. 15). 
When eqn.(78)-(82) are combined, the globally assembled finite-element solu- 
tion algorithm for the representative partial differential equation, system 
becomes 


i 
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The rank of the global equation system (83) is identical with the total number 
of node points on R 9R for which the dependent variable requires solution. 

5 Eqn. (83) is a first-order, ordinary differential system, and the matrix struc- 

j ture is sparse and banded. Solution of this system is obtained by COMOC using 

j| a predictor-corrector finite-difference numerical integration algorithm 
j: (ref . 14) . 

j A solution algorithm is required for the continuity equation, which is 

retained as eqn. (45) or (57) for boundary-layer or parabolic flows. When 
! retained, it describes an initial-value problem on pU 2 or p as a function 

t of x 2 , with x-j (x) and x 3 appearing as parameters. The solution approxi- 

! mation function need span only the transverse coordinate direction as 

i- ^ 

\ - 
f ’ 

| <5“ T (84) 

| . 

1 The matrix elements of Q are nodal values of pu| or p* *, their functional 

j; dependence requires solution of eqn. (45), (63), (66) or (65) along lines 

K (x.].x 3 ) or (x 2 ) equal a constant. Since each exists in standard form as an 
H ordinary differential equation, direct numerical integration or quadrature 
| i yields the required solution at node points of the discretization, 
j • 

I ; 

l 

i : 

J 

i 

i 
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COMOC COMPUTER PROGRAM 


The COMOC computer program system is being developed to transmit the 
rapid theoretical progress in finite element solution methodology into a 
viable numerical solution capability. In the course of establishing this 
general -purpose concept, several Variants o' COMOC have been developed for 
specific problem classes, including transient thermal and thermo-structural 
analysis, the two-dimensional transient Navier-Stokes equations, and the three- 
dimensional boundary-region, and parabolic Navier-Stokes equations. These 
initial forms have now been coalesced onto two advanced Variants. The Compu- 
tational Continuum Mechanics Variant is operational on IBM 360/ and 370/ com- 
puters; it serves as a research test bed and can solve multidisciplinary prob- 
lems characterized by a single differential equation description including 
transient and steady-state thermal analysis, thermo-el asto-statics, potential 
fluid flow, and magneto- and electro-statics. It contains automated data gen- 
eration features using curved iso-parametric two-dimensional elements to model 
non-regular shaped solution domains, and input/output graphics to facilitate 
data management and solution interpretation. The Navier-Stokes Variant is 
operational on the IBM 360/ and CDC 6000 series computers; it solves the 
multiple-differential equation descriptions characteristic of viscous fluid 
mechanics including transient-incompressible, and steady-compressible, 
multiple-species complete two-dimensional Navier-Stokes equations, as well as 
the two- and three-dimensional compressible, reacting, multiple-species 
boundary layer, boundary region and confined-flow parabolic Navier-Stokes 
equations. It also contains output graphics and an automated data generator 
for regular- shaped solution domains. An on-line restart feature allows the 
user to switch between bound ary- region and parabolic Navier-Stokes systems 
according to the requirements of the problem at hand. 

The finite element solution algorithm is utilized, as we observed in the 
previous section, to cast the original initial-valued, or pure-elliptic 
boundary-value problem description into large-order systems with a purely 
initial -value or algebraic character. COMOC then integrates or equation solves 
the discretized equivalent of the governing equation system. Initial distri- 
butions of all dependent variables may be appropriately specified or computed, 
and boundary constraints for each dependent variable can be specified on 
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arbitrarily disjoint segments of the solution domain closure. The solutions 
for each dependent variable, and all computed parameters, are established at 
node points lying on a specifiably nonregular computational lattice, formed by 
plane tri angulation of the elliptic portion of the solution domain. Each of 
the computational triangles is spanned by a linear approximation function used 
for all independent and dependent variables as well as all solution parameters. 

The COMOC system is being built upon the macros tructure illustrated in 
fig. 1. The main executive routine allocates core by means of a variable 
dimensioning scheme based upon the total degrees of freedom of the global prob- 
lem. The size of the largest problem that can be solved is thus limited only 
by the available core of the computer in use. The precise mix between depen- 
dent variables and parameters, and fineness of the discretization, is user- 
specifiable and widely variable. The input module serves its standard function 
for all arrays of dependent variables, parameters, and geometric coordinates. 
The discretization moduie forms the finite-element discretization of the ellip- 
tic solution domain and evaluates all required finite-element nonstandard 
matrices and standard-matrix multipliers. The initial ization module computes 
the remaining initial parametric data required to start the solution. The 
integration module constitutes the primary execution sequence of problem solu- 
tion, and primarily utilizes a highly stable, predictor- corrector integration 
algorithm for the column vector of unknowns of the solution. Calls to auxil- 
iary routines for parameter evaluation (viscosity, Prandtl number, source 
terms, combustion parameters, etc.) as specified functions of dependent and/or 
independent variables, as well as 'calls for equation solving algebraic systems, 
are governed by the integration module. The output module is similarly 
addressed from the integration sequence and serves its standard function via 
a highly automated array display algorithm. Both Variants of COMOC can execute 
distinct problems in sequence, and contain automatic restart capability to 
continue solutions. 













NUMERICAL SOLUTIONS IN AERODYNAMICS 


Both operational Variants of the COMOC system have been exercised to 
obtain finite element characterization of various aerodynamic flow field con- 
figurations. Linearized and complete potential flow solutions have been gener- 
ated by the Continuum Mechanics Variant, over curvilinear and airfoil configu- 
rations, to assess solution accuracy and to demonstrate features of finite 
element solutions. Viscous flow computations using the Navier-Stokes Variant 
include boundary layer prediction in a pressure gradient, parabolic Navier- 
Stokes solution of a trailing edge wake simulation, and a transient Navier- 
Stokes solution for decay of an aerodynamic shed vortex. 


Aerodynamic Potential Flow Solutions 

An informative problem geometry to evaluate finite element solution of 
linearized potential flow, eqn.(39), corresponds to subsonic, inviscid, iso- 
energetic flow over a wave-shaped wall, see fig. 2a. A linearized analytical 
solution can be obtained (cf. ref. 20, p. 458), for the ratio of wavelength 
(a) to amplitude (e) small, as 



for a sinusoidal wall described by 

y w (x) = E sin M ( 86 ) 

assuming a uniform inflow of at x = 0 and a frees tream Mach number of 
M m . An absolute accuracy assessment is possible, and the essential geometri- 
cal character requires a specified gradient boundary condition on a non- 
coordinate surface. The key computational output is surface velocity distri- 
bution, since from Bernoulli's law, eqn,{75), the corresponding incompressible 
surface pressure distribution can be determined as 

P(*,y w ) « P 0 - £ p„ u k u k 


(87) 


“ p o " 2 p ~ $, k -’k 
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The last form is obtained using eon. (28); formation of the indicated vector 
inner product, from the discretized computational solution for $ , is 
achieved using the finite element assembly operator ^ s see eqn.(81) and 
ref. 15. 

The existence of various symmetry planes in the problem domain can be 
effectively utilized to reduce the number of finite elements required to gen- 
erate a solution. Shown in fig. 2b-2c are economical domain closures and 
appropriate boundary condition specifications for a potential flow solution 
using either & or f . For the sample problem, 

U « 32 m/s 
00 9 

H = 0.2 

CO 

e/X = 0.025 

X = 2ir 

For the discrete approximation, the infinity boundary was chosen to lie at 
y = 60 e to insure that application of the zero gradient boundary condition 
was valid. The computational grid, consisting of 240 triangular finite ele- 
ments, was automatically generated from coordinate data describing vertex 
and mid-side nodal coordinates of two "super elements" as illustrated to the 
right in fig. 3a. A plot of the computed equi potential distribution appears 
in fig. 3b. 

Since the potential equation is elliptic, boundary condition specifica- 
tion is required on all domain boundaries. Indication of adequate discretiza- 
tion, to allow determination of accurate slope boundary condition representa- 
tion ,is noted both at the wall and at the upper boundary. The numerical and 
analytical solutions are compared in table 1, and the location of the maximum 
error is noted. To evaluate surface pressure, eqn.(87), the computed poten- 
tial function distribution is analytically differentiated using eqn.{80), and 
assembled over the entire domain to obtain scalar velocities and pressure 
along the wall. It is interesting to note that the level of velocity error 
is equal to that of the error in potential function, even though the former 
requires differentiation. The largest velocity error is 2.3%, indicating 
good accuracy for the coarse finite element mesh used. Since the pressure 
coefficient is related to the square of velocity, the corresponding maximum 
error in surface pressure is 5.3%. 
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b) Solution domain for 
scalar potential 


c) Solution domain for 
vector potential 


Figure 2 . Flow Over a Wave-Shaped Wall 






Transverse coordinate 



b) Equipotential distribution 
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Table 1. Potential Flow Over a Wave-Shaped Wall 


Potential Function - § 

Surface Velocity - 

r \ u k 

Analyti cal 

COMOC 

Error, % 

Analytical 

COMOC 

Error, % 

166.2 

166.7 

.3 

97.18 

97.15 

.03 

153.2 

153.8 

.4 

99.05 

98.16 

.91 

140.0 

140.7 

.5 

100.80 

99.54 

1.25 

126.5 

127.4 

.7 

102.40 

100.85 

1.51 

112.9 

113.9 

.9 

103.83 

102.12 

1.70 

99.1 

100.2 

1.0 

105.10 

104.33 

.74 

85.2 

86.4 

1.0 

106.19 

104.34 

1.75 

71.1 

72.3 

1.7 

107.11 

106.76 

.33 

57.0 

58.1 

1.9 

107.86 

108.16 

.28 

42,8 

43.7 

2.1* 

108,45 

109.42 

.90 

28.6 

29.2 

2.1* 

108.86 

110.54 

1.50 

14.3 

14.6 

2.1* 

109.11 

111.62 

2.30* 

0 . 

0 . 

0 . 

109.20 

110.55 

1.20 


^Maximum Local Error 


From the proven convergence character of the finite element solutions of 
linear (ref. 3) and non-linear (ref. 14,15) field problems, the S% error in 
pressure could be reduced to about 1% by a uniform doubling of the fineness of 
the employed discretization. Unfortunately, the computer CPU required to 
obtain the more accurate solution also increases dramatically (by a factor of 
up to 8). Bearing this in mind, computational experiments were conducted to 
ascertain the influence of particular selected discretizations and closure 
locations (especially the infinity. boundaries) on solution accuracy. For this 

i 

study, however, the full tensor potential flow equation was solved, eqn.(34), 
using a linear iteration algorithm and sequential update of the effective 
diffusion tensor, [6. . - c" 2 $,.] . The geometry selected is a symmetric 

t J * J 

NACA 0015 airfoil at zero angle of attack in subsonic flow in an inviscid wind 
tunnel, see fig. 4. 

Since the circulation is zero for this case, only half the airfoil geo- 
metry need be considered, and the Kutta condition is intrinsic. The flow 
domain was automatically discretized by COMOC into 192 triangular finite 
elements, from user specification of the nodal coordinates of a coarse dis- 
cretization consisting of three super elements, fig, 4 . The infinite boun- 
dary (tunnel wall) was set at about 13 airfoil thicknesses from the centerline. 
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FINITE ELEMENT 
DISCRETIZATION 


Figure 4 



SUPEk ELEMENT 
DISCRETIZATION 


Finite Element Discretization for NACA 0015 Airfoil 



i 


( 


Flow conditions correspond to air at - 0.1 and uniform velocity at the 
airfoil leading edge. Figure 5 presents the resultant computed pressure dis- 
tribution {0^) obtained from eqn. (87), along the airfoil surface compared 
with a conformal transformation solution (ref. 24, p. 324). While the solu- 
tion trends are consistent, the "peaky" finite element results indicate the 
selected discretization was inadequate for definition of the large velocity 
gradients occurring at the airfoil leading edge. This problem existed but 
was less apparent in the wavy -wall solution, due to the much lower leading 
edge curvature, see fig. 6. Even though both finite element solutions used 
approximately identical discretization fineness, redefinition and refinement 
is required for the airfoil leading edge. 

Several numerical experiments were conducted to evaluate discretization 
influences on solution accuracy. Equation (34) defines an elliptic boundary 
value problem in subsonic flow; it therefore is inappropriate to have either 
the airfoil leading or trailing edge coincide with the solution domain closure, 
even though the finite element algorithm solves for boundary nodes lying on 
symmetry planes (gradient boundary condition) along with interior nodes, see 
eqn. (83). While retaining the discretization fineness of fig. 4, additional 
regions extending one chord length upstream and downstream of the airfoil 
were added to the solution domain. The computational discretization, contain- 
ing approximately 480 finite elements, is shown in fig. 7, and was automatic- 
ally generated by C0M0C from the five super element description shown in the 
lower half of the figure. Upon viewing the results for this discretization, 
a non-uniform refinement of the computational zones upstream of the leading 
edge was defined for C0M0C (by two input number changes). This produced an 
approximate halving of the longitudinal span of the finite elements directly 
in front of the leading edge. The influence of these discretization changes 
on computed pressure coefficient distribution is shown in fig. 8. The addi- 
tion of the coarser discretization upstream produced even poorer agreement 
with the reference data and the unsatisfactory "peakiness" remains. Dramatic 
improvement is noted for the upstream refined grid solution, even though the 
discretization over the airfoil remains identical to that shown in fig. 4. 

Some additional refinement on the leading edge region would further improve 
agreement at a modest increase in computer cost. Well into the upstream 
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Figure 5. Computed Pressure Coefficient Distribution 
on NACA 0015 Airfoil. 
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region, note how the finite element solution anticipates the existence of the 
airfoil; this is the character of an elliptic boundary value problem. The 
refined grid solution is in excellent agreement with the reference values in 
the 20-95% chord region. However, note how the finite element pressure 
level returns to frees tream immediately downstream of the trailing edge, 
wherein the impact of the Kutta condition on the reference solution drives the 
coefficient to zero. Of course, in actuality, trailing edge regions exhibit 
significant viscous effects. These influences effectively blunt the trailing 
edge for a potential flow analysis, and a specific accounting of viscous/ 
inviscid interaction is required to accurately simulate the physics. 

Alternative forms of grid manipulation exist for improving solution accu- 
racy in leading and trailing edge regions. Global grid refinement, while 
effective, increases the order (size) of the solution matrix to be solved 
dramatically, thereby requiring more core storage. In addition, since solu- 
tion time varies as approximately the square of the number of nodes, use of 
very fine grids is economically impractical. Local refinement can retain a 
manageable number of elements by employing non-uniform distributions to place 
smaller elements in the high velocity gradient areas. This may become quite 
practical, but considerable numerical experimentation is required to optimize 
the procedure. For highly curved surfaces, like the leading edge region, 
higher-order curved finite elements could be used, but their efficiency for 
non-linear solutions remains to be quantified (see ref. 8, 9). For low sub- 
sonic flows where density is essentially constant, deVries and Norrie (ref. 5) 
suggested, but did not document an alternative approach involving global 
refinement on a local basis. This can be accomplished by using the orthogo- 
nality properties of the potential function and streamf unction to shrink the 
solution domain, i.e., translate the infinity closure nearer the airfoil. 

This method appears promising, since the total number of nodes in the solu- 
tion field remains small while establishing progressively finer grids near 
the airfoil. To evaluate this method, the previous case of low Mach number 
flow over a symmetric NACA 0015 airfoil at zero angle of attack was repeti- 
tively solved starting with the discretization shown in fig. 7. The distri- 
bution of streamf unction was first determined by solution of eqn.(40) in two- 
dimensional form. The results are shown in fig. 9. The ij; = ,15 streamline 
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Figure 9. Initial Potential Flow Computed Streamfunction Distribution. 
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Figure 10. Refined Grid Finite Element Solution Domain for 
NACA 0015 Airfoil. 




Figure 11. Refined Grid Potential Function Solution for 
NACA 0015 Airfoil. 
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Figure 12. Second Refined Finite Element Grid for NACA 0015 Airfoil 






was then chosen as the upper boundary for the next sequential evaluation; the 
solution domain was rediscretized as shown in fig. 10. This discretization 
employs the same number of triangular finite elements as did the original, 
but the average finite element span has been reduced by approximately six. 
Figure 11 shows the computed potential solution on the reduced domain, and 
fig. 12 illustrates the next reduced solution domain, . lerein the upstream 
and downstream boundaries lie on constant potential .evels as computed using 
the discretization of fig. 11. This procedure appears to provide a viable 
alternative for handling infinity boundary conditions associated with external 
incompressible flows. However, it remains to fully assess the accuracy im- 
provement that results in pressure coefficient distribution on the basis of 
cost effectiveness. 

Dealing with non-symmetric airfoils and/or angle of attack requires full 
discretization of multiply-connected solution domains. Practical solution of 
these cases with lift is of ultimate importance. The major difficulties asso- 
ciated with these solutions is again associated with the leading and trailing 
edges, as well as application of the Kutta condition. To evaluate the accuracy 
and stability of a finite element solution, the solution domain of fig. 7 was 
approximately doubled in extent. The finite element discretization, obtained 
using a 10 super element specification, is shown in fig. 13. The discretiza- 
tion is purposely non-symmetric, to evaluate solution accuracy; the resultant 
finite element mesh contains 355 elements and 213 nodes. A visually symmetric 
potential distribution was computed, as shown in fig. 14. Nevertheless, the 
discretization produced numerical solution differences for pressure coefficient 
on the upper and lower surfaces, as shown in fig. 15. The maximum differences 
are on the order of 15%; agreement with the reference data is essentially com- 
parable to that shown in fig. 8. A lift coefficient, , was computed to 
determine the error induced by these differences in pressure coefficient dis- 
tribution. The lift coefficient is computed by numerical integration of the 
upper and lower levels around the airfoil as (ref. 25). 

C L = t C P u " C P^ do: ( 88 ) 
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Figure 13. Complete Domain Finite Element Discretization 
for NACA 0015 Airfoil. 


Figure 14. Complete Computed Potential Function Distribution 
for NACA 0015 Airfoil. 
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Pressure coefficient 





The resulting determination for the NACA 0015 airfoil, at zero angle of 
attack, see Table 2, is -.00276, or an effective attack angle of -.02° . 

These results tend to indicate the acceptable accuracy in lift coefficient 
can be achieved using coarse and skewed discretizations. This is quite impor- 
tant, since for nonsymmetric shapes at angle of attack, it is essentially 
impossible to maintain a symmetric discretization. The magnitude of the 
error in C , hence C L , is directly associated with discretization fineness; 
as previously demonstrated, it can be controlled through grid refinement. 


Table 2. Pressure Coefficient Distribution on NACA 0015 
Symmetric Airfoil at Zero Angle of Attack, 
Non-Symmetric Finite Element Grid 


X 

C P U 

C Pt 

VS 

As 

i /V 

1.0 

.380 

.380 

0 . 

0 . 

0 . 

1.0333 

.035 

.180 

-.145 

.0333 

-.00240 

1.0667 

-.517 

-.298 

-.101 

.0333 

-.00410 

1.1 

-.442 

-.416 

-.030 

.0333 

-.00218 

1.1667 

-.473 

-.461 

-.012 

.0667 

-.00141 

1.2333 

-.491 

-.448 

-.043 

.0667 

-.00183 

1.3 

-.444 

-.425 

-.019 

.0667 

-.00210 

1.338 

-.440 

-.444 

.004 

.038 

-.00028 

1.4 

-.364 

-.368 

.004 

.062 

.00025 

1.488 

-.291 

-.292 

001 

.088 

.00022 

1.6 

-.210 

-.213 

.003 

.112 

.00022 

1.7 

-.128 

-.164 

- . .036 

.1 

.00200 

1.8 

-.033 

-.096 

.063 

.1 

.00495 

1.9 

.075 

.012 

.063 

.1 

.00630 

2.0 

.145 

.145 

0, 

.1 

.00320 


C L = - .00276 


As an example for an unsymmetric flow, the airfoil of fig. 13 was rotated 
to 15° angle of attack, and solution domain automatically discretized, as 
shown in fig. 16. The ability of the finite element procedure to use such 
non-uniform discretizations now becomes a distinct feature-, the capability to 
produce them automatically can reduce user input effort by several orders of 
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Figure 16. Complete Domain Finite Element Discretization for 
NACA 0015 Airfoil at 15° Angle of Attack. 
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magnitude. As embodied in COMOC, automatic refinement of minimal data input 
is accomplished over arbitrary geometric shapes. The algorithm involves quad- 
ratic functional representation over one, two and three-dimensional finite 
element spaces, using natural coordinate descriptions of the serendipidity 
finite element family (ref. 26). For the two-dimensional problems considered 
herein, the user-specified super elements are general triangular or quadri- 
lateral shapes, which may have curved sides not containing points of inflec- 
tion. Any two-dimensional domain or series of domains can be discretized by 
decomposition into quadrilateral and triangular shapes and specifying the 
number and type of elements to be generated in each. The level of user input 
effort can be even further reduced when geometric similarity -^f boundary 
regions exists for a class of problems. Airfoil analysis at various angles of 
attack in an unbounded stream exhibits such a similarity which can be exploit- 
ed. To illustrate the general case, consider the sequence in fig. 17. Uni- 
form inflow conditions are specified at the left end of the rectangular (or 
any other shape) box of fig. 17a. A vanishing normal gradient is appropriate 
along the top and bottom streamlines, while at outflow, proper surface orien- 
tation allows specification of constant potential. The airfoil shape, see 
fig. 17b, is input as specification of mean line coordinates and thicknesses. 
From these two specifications, super element coordinate data can be generated 
in a predetermined fashion for angle of attack, see fig. 17c-d, to serve as 
automatic discretizer input data. The addition of flaps or slats is concep- 
tually straightforward; options could be added to the super element genera- 
tion routine to allow their specification in terms of chord line coordinates, 
thickness distribution, and angle of attack. The combining of these tech- 
niques, coupled with automatic setting of the domain boundaries along lines 
of constant streamf unction and potential function, can provide a powerful, 
rapidly accessible and reliable tool for finite element analysis of general 
airfoil configurations. 

Aerodynamic Viscous Flow Solutions 

As indicated earlier, for conventional aerodynamic flows at low angle of 
attack, the inviscid flow field is separated from the airfoil by a generally 
thin region dominated by viscous effects. Therefore, the solution of the 
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a) Infinity Boundaries (f (t) ) b) Airfoil Description 
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c) Automatically Generated Initial 
Super Element Discretization 


d) Automatically Generated Streamline 
Super Element Discretization 


tn 

w Figure 17. Automated Finite Element Discretization for Potential 

Flow Analysis. 


two-dimensional boundary layer equations is an important requirement for com- 
putational simulation, and is readily accomplished within the developed 
Navier-Stokes Variant of COMOC. Since this computer program assumes all 
"parabolic" flows are three-dimensional, the dimensional -degeneracy of two- 
dimensional flow is obtained by employing a single column of finite elements 
spanning the boundary layer thickness, see fig. 18. The discretization extends 
into the freestream, where the inviscid flow is matched by a vanishing grad- 
ient boundary condition. The wall is assumed no-slip, and the lateral vanish- 
ing gradient boundary conditions yield the desired two-dimensional simulation. 
The character of the finite element solution of eqn. (45)-(49) can be evaluated 
for accuracy and convergence by comparison with solutions produced by finite- 
difference techniques and with a similarity solution for constant specific 
heat. The check case corresponds to a nominal Mach 5, laminar, two-dimensional, 
air boundary- layer flow over an adiabatic wall in a favorable pressure grad- 
ient. With the assumption of constant specific heat, the flow is isoenergetic 
and it is necessary only to solve the momentum equation and the continuity 
equation. The initial distribution for longitudinal velocity u^ is estab- 
lished from the similar solution for 3 = 0.5 and S = 0 (ref. 27). The 
initial distribution for u 2 is obtained iteratively, and Sutherland's law 
is employed to compute viscosity. 

The test case is initialized at x-j - 0.03 m downstream from the leading 
edge. The boundary-layer thickness at this station 6 Q is 0.0039 in , the 
local Mach number M e is 3.77 , the Reynolds number R e is 0,83 x 10 5 per 
meter, and the adiabatic wall temperature T w is 1000 K . Shown in fig. 19 
are the computed distributions of skin friction, local frea-stream Mach 
number, and boundary-layer thickness for the case of constant specific heat. 
These were obtained with two uniform finite-element discretizations corres- 
ponding to four and eight elements spanning the initial boundary- layer thick- 
ness, see fig. 18. The static pressure distribution P e U-j) is also presented 
for reference. The boundary-layer thickness has increased more than fourfold 
within the solution. Only small differences, on the order of about 2 percent, 
exist between the two solutions, the finer discretization producing a slightly 
larger skin friction and smaller local Mach number. Superimposed for compar- 
ison purposes are the results for the similar solution (ref. 27) and a 20-zone 
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Figure 18. Finite Element Discretization for Two-Dimensional 
Boundary Layer Flow. 
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Figure 19. Computed Supersonic Boundary Layer Parameters 
M = 5, Re.. = .83 (5) /m s 8 = 0.5. 



finite-difference solution obtained with the Von Mises coordinate transforma- 
tion. Agreement among the four solutions is excellent (within 2 percent) for 

skin friction. The similar solution for M lies between the COMOC and 

e 

finite-difference solutions, and overall agreement is within ± 3 percent. 

Shown in fig. 20 are computed velocity profiles at x-j/ ** 22.7 , which 
is about midway through the presented solution. For reference purposes, the 
initial longitudinal velocity profile is given with the node locations of the 
four-element discretization superimposed. Both finite element solutions pro- 
duce u-j distributions that are slightly more concave upward in the midregion 
in comparison with the similar or finite-difference solution. The eight-element 
COMOC solution lies closer to the similar solution in the region where the two 
finite element solutions differ. The finite-difference solution lies appre- 
ciably below both the COMOC and similar solutions near the free stream. The 
computed transverse velocities, which are also plotted in fig. 20, show only 
slight differences between the two discretization solutions. The trends of the 
COMOC solutions are in excellent agreement with the established precedures; 
unfortunately, since each method of solution is distinctly numerical, no abso- 
lute accuracy assessment is established. However, for an incompressible 
boundary- layer flow, absolute accuracy and convergence rates for the finite- 
element solution have been established to be close to theoretically predicted 
values (ref. 14, 28). 

As noted in the discussion of potential flow solutions on a NACA 0015 
airfoil, the computed pressure coefficient distribution (see fig. 8) indicates 
the (inviscid) surface velocity is retarded by an adverse pressure gradient on 
the interval 0.20 £ x/C £ 0.99 . Along the interval 0.85 < x/C £ 0.99 , the 
computed velocity is below the reference freestream value. The finite element 
solution computed a local acceleration, immediately downstream of the trail- 
ing edge, sufficient to return the surface flow to equilibration with the 
freestream. Of course, viscous effects would significantly modify these per- 
fect fluid results. The merging flow phenomena in the trailing edge region 
can be particularly complex, and a thorough computational analysis could be of 
considerable value, especially with regards to overall drag prediction. The 
finite element algorithm for solution of the boundary layer equations (45)-(49) 
merging into a parabolic Navier-Stokes equations (57)- (59) solution was 
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evaluated by COMOC. As the first step towards an airfoil trailing edge simula- 
tion, arbitrarily different, laminar boundary layer flows were allowed to 
develop on both sides of an infinitely thin flat plate. The plate was then 
computationally terminated, and the two boundary layer trailing edge velocity 
profiles (u-[ and u 2 ) were allowed to merge within a single, unbounded solu- 
tion domain. The test case corresponds to (essentially) constant density, 
isoenergetic flow of air at a nominal Mach number of 0.27 and a reference 
frees tream velocity, * 30 m/s . The finite element discretizations appear 
similar to that presented in fig. 18, (on both the top and bottom sides of the 
surface) and were of sufficient fineness to maintain ± 1 % solution accuracy on 
prediction of u-j . 

Initially, consideration was given to continued use of the boundary layer 
equation system downstream of the merging of identical upper and lower velocity 
profiles. Inconsistencies in the differential equation description were 
immediately encountered for u 2 , however. As described in the theoretical 
development section, the boundary layer equations cast solution for u 2 on an 
initial-value specification starting at the solid surface. This surface 
vanishes from the solution domain, immediately upon flow field merging down- 
stream of the plate termination. Since the flow ^rom the plate to the 
unbounded region is (assumed) smoothly continuous, in the latter the solution 
for u 2 must become a two-point boundary value problem (in two-dimensional 
space). This is possible only by retention of the x 2 momentum equation for 
solution of u 2 , and combining the continuity and vector momentum equation 
to yield an appropriately deterministic form for pressure prediction, see 
eqn. (62)-(66) . Computational experimentation using COMOC confirmed these 
theoretical musings. It was indeed impossible to obtain smooth transition 
from the boundary layer u 2 distributions while maintaining consistency of 
the frees tream values. 

These early results actually led to extension of COMOC to switch over to 
the appropriate parabolic Navier-Stokes system while undergoing a restart. 

For ensuing simulations, dissimilar upper and lower surface boundary layer 
profiles were allowed to develop from an initial slug (uniform) profile. 

Such a starting procedure minimizes input preparation; the boundary layer 
velocity profile that developed 1 0 2 m downstream of the simulated leading edge 
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agreed within ±1% of the Blasius solution at 1,0 in downstream at the same 
Reynolds and Mach numbers. Hence, arbitrarily dissimilar boundary layer pro- 
files were readily obtained to initiate the parabolic Navier-Stokes solutions. 
At restart, all program parameters used by COMOC to control integration were 
reset to their initial values. The two previously separate boundary layer 
solution domains were merged by adding the node at the plate into the inte- 
gration array, and specifying vanishing gradient boundary conditions for both 
u-j and u 2 everywhere about the closure of the newly defined domain. Initial 
evaluation for identical upper and lower initial velocity profiles showed that 
the COMOC- predicted results for u^ and u 2 remained exactly symmetric to 
five significant digits for merged solution domains up to 1.3 m long. The 
axial pressure gradient, computed using a pressure algorithm for eqn.(63), 
vanished to within 0.5% on U OT . The skew-symmetry on u 2 was exactly pre- 
served; the null value remained on the downstream projection of the plate 
throughout the solution. Subsequent evaluations utilized dissimilar upper and 
lower velocity profiles. The results of one such computation are shown in 
fig. 21. The plate terminus boundary layer initialization profiles are shown 
for x-j/<$ 0 = 267 ; other profiles are shown for various stations downstream of 
the trailing edge. Due to the original dissimilarity of the initial velocity 
profiles, note that the locus of the velocity profile minimum is concave 
upwards, yielding a modest overall curvature in the merged flow field. For 
this case also, the computed axial pressure gradient, eqn.{63), vanished to 
within 0,5% on U ra . The lateral pressure gradient was assumed to vanish 
identically. For non-flat plate flows, both axial and lateral pressure 
gradients can be induced by flow field curvature. Their computational pre- 
diction requires embodiment and check-out of solution algorithms for the 
various pressure descriptors, eqn. {65)-{66}. 

Computational predictions in laminar wakes are of limited interest in 
practical aerodynamics; the prime value of the discussed results is identi- 
fication of appropriate equation systems and proof of computational sta- 
bility. Even at low angle of attack, the initial laminar boundary layer 
flow will typically undergo transition and become fully turbulent before 
the trailing edge is reached. Hence, a far more practical computation for 
simulated airfoil wake flow would involve specification of a turbulent 
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Figure 21. Dissimilar Laminar Boundary Layer Merging Downstream 
of a Trailing Edge. 


closure relation, eqn. ( 5) » for solution of eqn«{45)-(48) and/or eqn. (57)- 
(58). Conventional turbulent boundary layers can be effectively predicted 
using a mixing length hypothesis (cf., ref. 29), wherein 

^12 = I U 1 ;2 I 

In eqn. (88), & is the mixing length conventionally defined as 

f k x 2 w 0 < x 2 < X6/k 

(A 5 iii X6/k < x 2 

where typically, k = 0.435 , X * 0.09 , & is the boundary layer thickness, 

and a) is the Van Driest damping coefficient, used to smoothly merge eqn. (88) 
into the sublayer molecular kinematic viscosity.. An alternative approach, of 
general applicability throughout multi -dimensional fluid mechanics, involves 
formation of the turbulent kinetic energy of the flow, and a closure hypothe- 
sis involving a suitable length scale (cf., ref. 18). Employing an order of 
magnitude analysis (cf., ref. 30), the turbulent kinetic energy is a point 
function which satisfies the same general partial differential system, eqn. 
(76)-(79), for which the finite element solution algorithm has been estab- 
lished. The required length scale can also be determined from solution of a 
differential equation (for example, from the dissipation function), or be 
hypothesized directly for geometrically simple flows. Results for finite 
element prediction of turbulent three-dimensional boundary region flows, 
using both closure techniques and-COMOC, are reported in refs. 31-32. These 
results indicate that extension of the finite element methodology to multi- 
dimensional turbulent flows, of impact in trailing edge wake flow predictions, 
will be directly accomplished. 

The final aerodynamic problem class, for which the finite element solu- 
tion procedure has been evaluated, corresponds to highly rotational vortex 
flows. Two problem classes of immediate potential applicability include, 

1) analysis of the persistence of trailing vortex streets as generated by the 
lift distribution along highly loaded airfoil configurations typical of 
current wi de-body jets, and 2) prediction of the penetration and turbulent 
decay of high energy jet flows injected non-tangentially into a crossflow as 
might occur in thrust vector control for VTOL aircraft. In each case, the 
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vortex structure will decay as a function of time and the local effective 
diffusion and convection effects. The vortex center will move as a function 
of cross flow velocity, initial strength and direction, and the magnitude of 
the tangential velocity at the juncture of contra-rotating vortex pairs. A 
complete numerical simulation typically requires solution of the full form of 
the Navier-Stokes equations. A cursory evaluation of the finite element solu- 
tion procedure has been explored for the second problem. An initially circu- 
lar jet, injected subsonically and non-parallel into a cross wind, is struc- 
tured as a contra-rotating vortex pair (cf. , ref. 33), For sufficient dynam- 
ic pressure, the jet penetrates into and eventually turns parallel to the 
impressed crossflow. The jet vortex structure provides the “elastic stiffness" 
essential for the impressed flow to pass around the initial jet as if it was a 
aerodynamic stream tube. As the jet turns and travels downstream, the impressed 
crossflow influences the dissipation of the initial highly-rotational structure 
by convection and diffusional processes,. 

Determination of a streamtube structure, as a function of initial vortex 
strength and assumed crossflow, was made using COMOC, by modeling the three- 
dimensional problem as a transient two-dimensional configuration described by 
the Navier-Stokes system written on stream function and vorticity, eqn.(69)-(70). 
A sketch of the two-dimensional vortex structure is shown in fig. 22a; the 
indicated counter-rotation deflects an onset flow, directed anti -parallel to the 
X 2 coordinate direction, around the centroidal region. The smoothed distri- 
bution of vorticity along the horizontal symmetry plane is shown in fig. 22b. 
For the computational simulations this was replaced by point sources of vor- 
ticity set equal to ±w appropriately located at x-j = ±x Q . Existence of 
the mirror symmetry plane, x-j = 0 , combined with use of a "sufficiently 
large" solution domain allows specification of vanishing vorticity everywhere 
on the selected domain closure for solution of eqn.(70). The mirror plane is 
also a streamline, as is the right side boundary in fig. 22a. The value of 
streamf unction taken on the latter is determined by the specified magnitude 
of the impressed crossflow, assumed anti-parallel to the X 2 axis at inflow 
(only). The resultant solution of eqn.(69) for streamf unction determines the 
actual inflow and outflow velocity distributions with use of the vanishing 
gradient (parallel flow) boundary condition for streamf unction, 
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a 0) = a ( 3 ) = o , eqn = (78). With these imbedded symmetry properties, only 
one-half the domain shown in fig. j 22a is required for determining a solution; 
the right half was discretized usjing 200 triangular-shaped finite elements and 
121 nodes. The cavity was presumed filled with air initially at rest, and the 
single vortex pair, of variable strength, was located at x-j = ±20% h , where 
h is the half-width domain span, fig. 22a. 

Shown in fig. 23 are the COMOC computed distributions of inital mass flow 
within and through the solution domain for three specified values of initial 
vorticity, to 0 . The displacement effect of the point source is graphically 
apparent. For zero specified crossflow, the computed initial streamf unction 
forms closed contours, as shown in fig. 24 for different initial locations for 
the vortex pair. (The spurious wiggles result from use of high order spline 
fits to sparse data by the plot package.) A continued solution of eqn.(69)- 
(70), to predict the decay of the initial point vortex pair, requires speci- 
fication of an appropriate turbulence closure relation, eqn.(5). An accurate 
simulation would also require that the impressed crossflow be time varying. 
This would be readily accomplished in COMOC by altering the value of stream- 
function on the right closure segment, fig. 23, as a specified function of 
time. 
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Figure 23. Initial Streamfuncti on Distributions for Cross Flow Onto a Point Vortex Pair. 





CONCLUSIONS 


The theoretical and computational results of this study, on application 
of finite element solution methodology to configuration analysis in low speed 
aerodynamics, indicate significant potential for the procedure. The approach 
taken has established appropriately complete (i.e., non-linear) differential 
equation descriptions for the several distinct fluid flow problem classes of 
impact in aerodynamics. The developed finite element algorithm is universally 
applicable to each description. The developing CQMOC computer program, which 
embodies this algorithm, has verified the overall concept of a powerful, ver- 
satile general-purpose code for computational low speed aerodynamics. The 
generated numerical solutions in inviscid potential flow have introduced and 
evaluated various techniques for error control within the constraints of 
computational and input- requirement economics. Diverse viscous flow predic- 
tions, including a transition of differential equation system during a problem 
solution, indicate broad-base applicability of the algorithm and its computa- 
tional embodiment. Continued progress on development of appropriate pressure 
solution algorithms for parabolic flows, coupling of the inviscid and viscous 
solutions, and expanded input/output graphics should yield a highly useful, 
versatile and user-oriented design tool to supplement and supplant detailed 
wind tunnel evaluation of complex low-speed aerodynamic systems. 
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